Long-time saturation of the Loschmidt echo in quantum chaotic billiards 
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The Loschmidt echo (LE) (or fidelity) quantifies the sensitivity of the time evolution of a quantum 
system with respect to a perturbation of the Hamiltonian. In a typical chaotic system the LE has 
been previously argued to exhibit a long-time saturation at a value inversely proportional to the 
effective size of the Hilbert space of the system. However, until now no quantitative results have 
been known and, in particular, no explicit expression for the proportionality constant has been 
proposed. In this paper we perform a quantitative analysis of the phenomenon of the LE saturation 
and provide the analytical expression for its long-time saturation value for a semiclassical particle 
in a two-dimensional chaotic billiard. We further perform extensive (fully quantum mechanical) 
numerical calculations of the LE saturation value and find the numerical results to support the 
semiclassical theory. 



I. INTRODUCTION 

In his seminal 1984 paper [l(, Peres studied the stabil- 
ity of motion of quantum systems with respect to small 
perturbations of the Hamiltonian. He discovered that 
the quantum motion of a system, whose underlying clas- 
sical dynamics is chaotic, is more unstable than that of a 
system, whose dynamics is regular in the classical limit. 
The quantity introduced by Peres, presently known as 
the Loschmidt echo (LE) or fidelity, has been the subject 
of thorough theoretical and experimental research in the 
fields of quantum chaos and quantum information 0, Q • 

The LE, defined as 
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quantifies the "distance" (in the Hilbert space) between 
the state e~ lHt / h \4>o) , resulting form the initial state |<^o) 
in the course of evolution through a time t under the 

Hamiltonian H, and the state e~* Ht / h \4>o) obtained by 
evolving the same initial state through the same time 
t, but under a slightly different, perturbed Hamiltonian 
H . The LE, by construction, equals unity for t = and 
typically decays further in time. A variety of different 
decay regimes - the most prominent ones being the Lya- 
punov [J, Fermi- Golden- Rule [H, [H, and the perturbative 
@j B regime - have been found in chaotic systems with 
various Hamiltonians and Hamiltonian perturbations. In 
this paper, however, we address the property of the LE 
generally shared by all (Hermitian) chaotic systems: the 
saturation of the decay at long times. 

Peres provided in his original work [lj a qualitative 
(order-of-magnitude) estimate for the value M m of the 
LE saturation in chaotic systems. He argued that for 
small enough perturbations 



where N is the number of eigenstates (of the unperturbed 
Hamiltonian H) that are significantly represented in the 
initial state \4>o)- In other words, N is the size of the 
effective Hilbert space that is required for a reasonable 
description of the time evolution of the initial state. 

The phenomenon of the LE saturation has been pre- 
viously addressed in the literature from numerical Q 
and analytical [9] perspective, and the validity of the 
Peres' argument, Eq. ([3]), has been verified. However, 
no explicit expression for the proportionality constant in 
Eq. ([3]) has been proposed. Our work complements the 
theory of the LE in chaotic systems by providing the 
(previously missing) proportionality constant. 

In this paper we present the semiclassical analysis of 
the LE at long times, and derive an expression for the 
LE saturation value. Our result, while in agreement with 
Eq. ([3]), constitutes a quantitative estimate of M^. The 
system treated in this paper is a two-dimensional, quan- 
tum billiard that exhibits chaotic dynamics in the clas- 
sical limit. The key method underlying our analytical 
calculation, however, is not restricted to billiards and 
can be generalized to a wider range of chaotic systems. 
We further perform numerical simulations of the time 
evolution of an initially localized, Gaussian wave packet 
in a chaotic billiard, and compute the saturation of the 
LE due to a perturbation caused by a deformation of 
the boundary. The results of the numerical simulation 
strongly support our analytical predictions. Finally, we 
conclude the paper with a discussion and final remarks. 



II. SEMICLASSICAL APPROACH 

We consider the time evolution of a quantum particle 
moving inside a two-dimensional ballistic cavity - a quan- 
tum billiard. In this paper we only consider hard-wall 
billiards whose underlying classical dynamics is fully hy- 
perbolic [l0( • The initial state of the particle is assumed 
to be the coherent state 
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with a quantifying the dispersion of the wave packet, 
and r and po representing respectively the average po- 
sition and momentum of the particle. The dispersion a 
is assumed to be small compared to the linear size of 
the billiard for the wave function to be normalizable to 
unity. We further define the de Broglie wavelength of the 
particle as 
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Po 
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where po = |Po| is the magnitude of the particle's momen- 
tum. (Hereinafter we denote the magnitude of a vector 
by its corresponding symbol in italics.) 

The time evolution of the initial state in the un- 
perturbed system with the Hamiltonian H is given by 
</>t(r) = J dr' Kt(r,r')(f>o(r'). In the (short-wavelength) 
semiclassical approximation the propagator Kt(r, r'), for 
a two-dimensional system, can be written as [111 ] 
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Here, Sy denotes the action integral along the classi- 
cal path 7' leading from the position r' to r in time 

t, and D Y = |det(-<9 2 SV/<9r<9r')| 1/2 e 
Maslov index /.ty 
Ref. 0) 



W 2 with the 
Then, in the limit (see Appendix A of 
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with It, being the Lyapunov length of the billiard, the 
action integral Sy can be linearized about the trajectory 
7(ro — > r, t) connecting the wave packet center ro and the 
point r in time t: Sy ~ / 



-n (i) 
p 7 



(r' — r ), where p 7 is 



the initial momentum on the trajectory 7. Using this ac- 
tion linearization and performing a Gaussian integration 
over the initial point r' one obtains the semiclassical ex- 
pression for the time-dependent wave function evolving 
under H 0: 
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The wave function <fit(j) corresponding to the time evo- 
lution under the Hamiltonian H of the perturbed system 
is given by an equation analogous to Eq. (JSJ) with the tra- 
jectories 7(r — > r, f) replaced by 7(r — > r, t) satisfying 
the classical evolution corresponding to H. 
The LE amplitude, 0(t) = (<fit\(j)t), is given by 
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The expression for the LE, then, being the product 
0*(i)0(t), with the asterisk denoting the complex conju- 
gation, involves two integrals over the final points, say r 
and f , over four sums over trajectories, two corresponding 
to the perturbed system and two to the unperturbed one. 
The integrand, in general, is a rapidly oscillating func- 
tion of r and f ; therefore, only the trajectories with the 
overall phase difference smaller than h give a finite con- 
tribution to the integral. Considering trajectories such 
that S-y Sty leads to exponentially decaying regimes 
of the LE Therefore, non-decaying contributions to 
the LE (responsible for the LE saturation) can only re- 
sult from trajectories that are close in action and belong 
to the same Hamiltonian. This imposes a restriction on 
the possible configurations of the trajectories of interest, 
namely r w f. This makes it convenient to make the 
following transformation to the new integration coordi- 
nates: Q = (r + f)/2 and q = r f . Then, following the 
procedure above, we linearize the four trajectories enter- 
ing the expression for the LE about the same final point 
Q to obtain 
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where AS = (5 7 — 5 7 
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Py ) ' Q/2) ah the four paths (7, 7', 7 and 7') connect 
r and Q in time t, with two of them (7 and 7') cor- 
responding to the unperturbed Hamiltonian H and the 
other two (7 and 7') to the perturbed Hamiltonian H; 
here p( f ) denotes the final momentum (at the end point 
Q) on the corresponding classical path. The integrand 
in Eq. ([TO]) is still a rapidly oscillating function of Q - 
AS is generally much greater than h and is sensitive to 
Q - unless the paths 7, 7' , 7 and 7' are correlated. In 
the diagonal approximation [131 ] the main contribution to 
the Q-integral comes from such terms in the sum that 
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One group of such terms, de- 



fined by the identification 7 = 7 and 7' =7', is responsi- 
ble for the (generally exponential) time-decay of the LE 
[H, and leads to a vanishing contribution at long times. 
The other group, defined by the identification 7 = 7' and 
7 = 7', gives rise to a term surviving in the limit t — > 00 
and, therefore, provides the leading order contribution to 
the LE saturation value. Thus, identifying the trajecto- 
ries of the unperturbed (7 = 7') and perturbed (7 = 7') 
Hamiltonian we obtain 
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In order to evaluate the double sum in the right hand 
side of Eq. (TiTj) we utilize the sum rule [Tjj] 

J2 |i? 7 | 2 /(r,P«;r',p( f )) 

7(r— >r',i) 

= J dp J dp'^(r,p;r',p')/(r,p;r',p'),(12) 

where Vt(r, p; r', p') = <S(r 7 (t) — r')d(p 7 (t) — p') is the 
classical phase-space probability density for a trajectory 
7 = {(r 7 (r),p 7 (r)),T G [0, t]} starting from the phase- 
space point (r 7 (0), p 7 (0)) = (r, p) to end at the point 
(r 7 (t), p 7 (i)) = (r', p') while evolving under the Hamilto- 
nian H through time t. Then, since dealing with chaotic 
Hamiltonians and long times, we replace the probability 
distribution Vt by its phase-space average 

6(H(r',p')-H(r,p)) 



V(v,p;r',p') 



Q(H(T,p)) 
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where Cl(E) is the phase-space volume of the energy shell 
H(r,p) — E. For the case of two-dimensional billiards 
Q,(E) = 2irmA, with m being the mass of the particle 
and A the billiard area, so that in view of Eqs. (| 1 2|) and 
(fT3| the long time (t — > oo) limit of Eq. (fTTjl reads 



dQ / dq 



dpdp'dpdp' S(H{r ,p) - H(Q, p')) 



xS(H(v ,p) - H(Q, p')) exp |-(p' - p') • q 
-^[(p-Po) 2 + ( P -Po) 2 ]}. 
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We now assume that both Hamiltonians can be writ- 
ten as p 2 /2m + V(r) and perform the integration in the 
right hand side of Eq. (Tj"4"|) as follows. The q-integration, 
with the integration limits extended to R 2 , results in 
{2ttH) 2 8(p' — p'). Consequently integrating over p' and 
p' we obtain 
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x<5(£(r ) - E(Q) + (p 2 /2m - p ,2 /2m)) , (15) 

where £(r) = V(r) — V'(r). Now we assume that the 
perturbation is small compared to the kinetic part of the 
Hamiltonian. Alternatively, one may consider perturba- 
tions of the Hamiltonian produced by deformations of the 
billiard boundary [l2|, [l5j]; it is the perturbation of the 
latter type that we use in our numerical experiments of 
the following section. Thus, assuming £ = 0, we have 



exp (~ [(p - p ) 2 + (p - Po) 2 ] | 
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Now, we use the integral representation of the 5- 
function, 5(p 2 — p 2 ) = (27r) _1 J d£exp(i£p 2 — it;p 2 ), and 
perform the Gaussian integration over p and p (eventu- 
ally doing the variable change x = (,h 2 /a 2 ) to get 
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where a = (poa/h) 2 = (27rcr/A) 2 , and Iq is the zeroth 
order modified Bessel function of the first kind. In the 
limit a 1 (or A -C u), which is in agreement with 
Eq. 0, the asymptotic form Io(a) w (2na)~ 1 ^ 2 exp(o) 
yields 
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Equation (fT8|) constitutes the central analytical result of 
our paper. 

It is easy to see that the original argument by Peres, 
see Eq. (|3]), is in perfect agreement with Eq. (| 18|) derived 
in the semiclassical approximation. Indeed, the number 
of Hamiltonian eigenstate required to properly describe 
the time evolution of the initial wave packet, given by 
Eq. (g]), can be evaluated as N = Vl(E)AE/(2kK) 2 . Here, 
as above, Q,(E) — 2nmA is the phase-space volume of 
the energy shell at the average energy E = p\j2m of the 
particle, and AE = poAp/m is the energy dispersion of 
the initial state. Estimating the momentum dispersion as 
Ap « 2V2H/ a we obtain the following expression for the 
number of the eigenstates: N w 2\[2Aj\a (and therefore 
Afoo w 2n~ 1 / 2 N~ 1 ). In fact, due to certain arbitrariness 
in determination of Ap the size of the effective Hilbert 
space N is not properly defined. This difficulty points 
to a drawback of the original formulation of Eq. ([3]). On 
the contrary, Eq. (TT5)) gives the LE saturation value in 
terms of well defined system parameters, A, cr, and A, 
and, therefore, provides a quantitative estimate for M^. 

In the following section we demonstrate that the semi- 
classical predictions of Eq. (TT8"]) are in agreement with 
the time saturation of the LE observed in numerical ex- 
periments. 



III. NUMERICAL SIMULATIONS 

In order to support our semiclassical calculations we 
have performed numerical simulations of a quantum par- 
ticle moving inside a desymmetrized diamond billiard 
(DDB). The DDB is defined as a fundamental domain of 
the area confined by four intersecting disks centered at 
the vertices of a square. The billiard is fully chaotic [l6[ 
and has been previously considered for studying various 
aspects of quantum chaos [l2l . ITU IT7T |. In our numer- 
ical experiments we used the piston-like boundary de- 
formation [T2I ] as the perturbation of the Hamiltonian. 
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The numerical method that we used for propagating the 
particle's wave function in time is the Trotter-Suzuki al- 
gorithm [Isj l; Reference [l2j provides further details on 
the billiard systems, Hamiltonian perturbation and wave 
function time propagation. 
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FIG. 1: (Color online) Time decay of the Loschmidt echo 
in the desymmetrized diamond billiard with a boundary de- 
formation for initial wave packets of de Broglie wavelength 
X — 4ir and dispersion a — 9. Time is given in units of 
the free flight time tt of the corresponding classical particle. 
The thin (blue) line shows an individual LE decay curve re- 
sulted from a single numerical experiment. The thick (red) 
line represents the result of an averaging over 3 individual 
decay curves obtained for different positions ro of the initial 
wave packet. 

In our simulations the initial state of the quantum par- 
ticle is given by Eq. ((4} . The blue line in Fig. Q] shows a 
typical LE decay curve obtained in an individual numeri- 
cal experiment with the initial wave packet of the disper- 
sion a — 9 and de Broglie wavelength A = 47r; the area 
of the billiard A « 1.51 x 10 5 . Time is given in units of 
the free flight time U of the counterpart classical billiard, 
i.e. t/tf is the number of bounces of the corresponding 
classical particle. The red line in Fig. [1] is the result of 
the averaging of the LE over 3 individual decay curves, 
each of which was obtained by propagating a wave packet 
centered about a different spatial point ro inside the bil- 
liard domain. (The wave packet centers were chosen such 
that the three initial states had negligible overlap with 
one another). Those were the average LE decay curves 
that we used to determine the LE saturation value and 
standard deviation - red dots and error bars in Fig. [2]- 
for the initial quantum state with particular values of the 
dispersion and de Broglie wavelength. 

In Fig. [2] we compare the semiclassical estimate for 
Mqo, given by Eq. (fT8|) . to the saturation values obtained 
from the numerical simulations. The top (bottom) figure 
shows the dependence of M M on the dispersion a (de 
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FIG. 2: (Color online) Top figure: The Loschmidt echo sat- 
uration value Moo as a function of the dispersion a of the 
initial wave packet for a fixed de Broglie wavelength A = 4ir. 
Bottom figure: Moo as a function of A for a — 9. In both fig- 
ures the billiard area A ps 1.51 x 10 , and the blue dashed line 
represents the LE saturation value as predicted by Eq. (JTSJ). 



Broglie wavelength A) for A = 4-7T (a — 9) in the billiard of 
the area A w 1.51 x 10 5 . The red dots together with the 
error bars represent the numerically observed values of 
Mqo ; the blue dashed lines are plotted in accordance with 
Eq. (fTg|) . We stress here that no free (fitting) parameters 
have been used in producing the theoretical lines: the 
slopes of the lines are entirely fixed by Eq. (fT8]) . 

Finally, to give an idea of the scale of the numerical 
simulations of this section we note that obtaining an indi- 
vidual LE decay curve, such as the blue curve in Fig. [TJ 
requires more than 8 days of computational time on a 
high-end (2.8GHz, 2GB RAM) computer. Each data 
point in Fig. [5] is a result of the averaging over 3 such 
individual decay curves. Therefore, 39 individual decay 
curves were obtained to produce the numerical data pre- 
sented in Fig. [2l amounting to approximately 312 days 
of (single-processor) computational time. 



IV. CONCLUDING REMARKS 

In this paper we have used the methods of the semi- 
classical theory to derive an explicit expression for the 
value of the long-time saturation of the LE, M^, in 
two-dimensional chaotic billiards. Our quantitative re- 
sult agrees with the early qualitative argument [l[ that 
the LE saturates at a value inversely proportional to the 
effective size of the Hilbert space of the system; our cal- 
culation provides the previously missing proportionality 
factor. 

In order to support our analytical predictions we have 
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performed careful numerical simulations of a quantum 
particle moving in a chaotic billiard. In these simula- 
tions a deformation of the billiard boundary played the 
role a Hamiltonian perturbation. The decay of the LE 
was observed until times long enough to reliably deter- 
mine Moo, and a proper ensemble averaging (over the 
initial position of the quantum particle) was performed 
to improve the accuracy. The numerically obtained val- 
ues of the LE saturation were found in a good agreement 
with the theory. 

The central aspect of our semiclassical calculation is 
the pairing (in the sense of the diagonal approximation) 
of trajectories that belong to the same (perturbed or 
unperturbed) Hamiltonian. Those are these trajectory 
pairs that render the time-independent contribution to 
the LE in addition to other, exponentially decaying con- 
tributions resulting from different trajectory pairs. Here 
we note that the pairing of trajectories considered in this 
work has been previously studied in the context of the 
fidelity fluctuations @ and the survival probability decay 
in open chaotic systems [l9j . 

We also note that although the phenomenon of the 
long-time saturation of the LE has been previously dis- 
cussed in the literature [8, 9] it has never been subject to 
a thorough analytical and/or numerical study. In partic- 
ular jthe numerical simulations of the quantum Lorentz 
gas [8( correctly demonstrated the inverse proportional- 
ity of Mqo to the billiard area, A, while misleadingly sug- 



gesting its linear dependence on the square of the wave 
packet dispersion, a 2 , along with independence of the de 
Broglie wavelength A. Reference [H, on the other hand, 
correctly outlined the semiclassical derivation of the di- 
rect proportionality of to the effective Plank con- 
stant, but did not present an explicit form of the propor- 
tionality coefficient. Thus the present paper bridges the 
gap by providing a quantitative analytical expression for 
the LE saturation value and, consequently, verifying the 
expression by means of extensive numerical simulations. 

As the final remark we would like to point out that 
the present semiclassical approach to the phenomenon 
of the LE saturation is only valid in the long time limit 
and in the regime of weak Hamiltonian perturbations. In 
general, however, the LE saturation value will depend on 
a (properly defined) perturbation strength. It is not yet 
clear to us how this dependence can be described by the 
semiclassical theory. 
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